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1.  INTRODUCTION. 
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Suppose  that  two  (possibly  dependent)  point  processes  are  observed 
simultaneously  over  a  period  of  time,  yielding  observations  at  times 

<  •••  <  for  the  first  process,  and  at  times  <  I$2  <  •••  <  B^ 

A  B 

for  the  second.  Such  data  arises  in  many  contexts,  and  it  is  often  of 

interest  to  discover  and  quantify  the  association  between  the  two  processes. 
Two  fields  in  which  this  situation  occurs  are  neurophysiology  and  reli¬ 
ability  theory,  from  which  the  following  four  examples  are  drawn. 

In  neurophysiology,  point  processes  arise  as  the  impulse  times  of 
neurons.  An  example  of  a  hypothesis  testing  problem  arises  when  one  wants 
to  determine  whether  or  not  the  impulse  times  of  two  neurons,  say  A  and  B, 
are  associated.  Of  course,  if  an  association  is  discovered,  it  is  usually  of 
interest  to  identify  the  nature  of  the  relationship  between  the  two  neurons; 
e.g.  is  A  driving  (inhibiting)  B?  If  so,  by  how  much?  This  is  typically  a 
harder  problem. 

A  second  example,  one  in  which  estimation  is  necessary,  is  given  by  the 
following  situation.  An  animal  is  to  be  taught  (or  trained)  to  perform  a 
certain  task.  Now  consider  two  connected  neurons  which  are  essential  in  the 
performance  of  this  task,  and  simultaneously  record  their  firing  patterns. 
Before  the  learning  has  taken  place,  the  neurons  will  in  general  affect  each 
other.  The  problem  is  to  determine  the  way  in  which  this  dependence  is 
altered  after  the  learning  has  taken  place. 
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The  next  two  examples  are  taken  from  reliability  theory.  Consider  a 
system  which  contains  two  components,  A  and  B,  and  for  which  A  is 
replaced  according  to  an  age  replacement  policy:  the  unit  is  replaced  at 

.  be  the 

successive  times  at  which  the  component  is  replaced  without  having  failed, 


failure  or  at  time  T,  whichever  comes  first.  Let  Ai  <  ^2  < 


and  let  B^  <  B2  <  • • •  be  the  times  of  failure  of  component  B.  It  may  be 
of  interest  to  determine  if  the  point  processes  1A^}  and  {B^}  are  depen¬ 
dent,  and  in  particular  to  determine  if  the  age  replacement  policy  is  bene¬ 
ficial  (or  harmful)  to  component  B. 

As  a  second  example,  consider  a  series  system  of  two  subsystems  A 
and  B,  where  A  and  B  are  parallel  structures  of  k^  and  kg  components, 
respectively  (kA  >_  2,  kfi  2).  Assume  that  any  failed  component  is  repaired, 
and  that  repair  time  is  negligible.  Thus,  if  any  component  fails  the  system 
continues  to  operate.  Let  A^  <  <  •••  and  B^  <  B^  <  •••  denote  the 

times  at  which  failures  occur  in  subsystems  A  and  B,  respectively. 

Then,  these  form  a  bivariate  point  process,  and  it  may  be  of  interest  to 
determine  if  the  failure  times  of  subsystems  A  and  B  are  dependent. 

In  this  article,  we  describe  and  discuss  certain  graphs,  plots,  as 
well  as  more  formal  methods  that  can  assess  the  dependence  between  point 
processes.  Specifically,  these  methods  indicate  whether  or  not  the  likeli¬ 
hood  of  an  A-point  is  increased  (decreased)  after  the  occurrence  of  a 
B-point.  The  techniques  are  illustrated  on  simulated  data.  Although 
bivariate  point  processes  arise  in  many  fields,  we  emphasize  the  applica¬ 
tions  in  neurophysiology. 

We  hope  to  illustrate  the  techniques  described  here  on  data  drawn  from 
real  neurons  later. 
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2.  THE  FUNCTION  XA|fi(*)  AND  THE  CROSS-CORRELATION  HISTOGRAM. 

We  proceed  at  a  heuristic  level.  A  rigorous  approach  and  formal 
definitions  of  the  terms  used  below  may  be  found  in  Daley  and  Vere-Jones 
(1972).  The  book  by  Cox  and  Isham  (1980)  is  a  good  guide  to  the  litera¬ 
ture  and  contains  all  the  necessary  information  to  formalize  what  we  do 
here. 

We  assume  that  the  points  {A^}  and  {B  }  form  a  bivariate  point 
process  that  is  stationary  and  orderly.  Intuitively,  the  stationarity 
assumption  is  that  the  process  has  been  going  on  for  a  long  time  and  is 
in  steady  state.  The  orderliness  condition  is  that  each  univariate  process 
has  no  multiple  occurrences.  For  a  Borel  subset  of  the  real  line  S,  let 
Na(S)  and  Ng(S)  denote  the  number  of  points  A^  and  B ^ ,  respectively, 
that  lie  in  S.  We  will  use  N^(s,t)  to  denote  ^((s.t]),  for 
i  =  A,B.  We  may  view  NA(*)  and  Ng(*)  as  random  Borel  measures. 

We  define  the  rates  ><A  and  Ag  as  follows: 

(1)  XA  =  lim  £  P{N  (t,t+h)  >  0}  , 

h-*0 

with  a  similar  definiton  for  Ag.  By  the  stationarity  assumption,  A 
and  Ag  are  independent  of  t.  The  existence  of  the  limit  in  (1)  was 
first  proved  by  Khintchine  (1960).  Korolyuk's  Theorem  (see  Leadbetter, 
1968)  states  that  if  t^  <  then 

(2)  E  NA(trt2)  «=  VW  • 


3 


Similarly  for  the  B  process.  A  simple  consequence  of  the  ergodic  theorem 
is  that  with  probability  one, 


N  (0,T) 

(3)  -  -  X.  as  T  , 

T  1 

for  i  =  A,B.  This  gives  a  third  way  of  thinking  about  the  rates  X * 
and  XB. 

The  discussion  so  far  relates  only  to  the  processes  %(•)  and  Ng(*), 
taken  one  at  a  time.  To  see  how  the  processes  affect  each  other,  we  define 
the  following  quantities: 

(4)  Xa!r(u)  =  lim  77  P^N. (t+u, t+u+h) | N- ({ t } )=1 }  for  -«  <  u  <  »  . 

A|B  h-K)  h  A  B 

Thus,  roughly  speaking,  Sives  the  infinitesimal  probability  of  an 

A-point  u  units  after  a  B-point.  The  stationarity  assumption  implies  that 
Xa]b(u)  i-s  independent  of  t.  We  may  also  define  the  function  X  g  |  ^  ( • )  * 
but  it  is  simple  to  see  that 


X 


B 


XA  AB|A^"U)  f°r  <  U  <  °°  » 


so  that  it  suffices  to  consider  \|b^u^’  3S  as  we  consi^er  both 

positive  and  negative  values  of  u. 

The  function  \]g(*)  has  been  considered  in  different  forms  and 
contexts  by  Cox  (1965),  Cox  and  Lewis  (1972),  Brillinger  (1976)  and 
Griffith  and  Horn  (1963),  among  others.  If  the  processes  N^(#)  and 
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Ng(*)  are  independent,  we  clearly  have 


Xa|B(u)  =  XA  for  a11  u  * 

so  that  X^jg(*)  can  indicate  deviations  from  independence.  The  following 
two  examples  of  very  simple  neuronal  networks  serve  to  illustrate  this 
point . 

Example  1:  Two  independent  neurons. 

N.  (• )  and  Ng(*)  are  independent.  Then,  as  was  just  mentioned, 

(5)  >A|B(u)  =  for  a11  n  ’ 

so  that  in  particular,  X^C*  )  is  constant. 

Example  2:  A  network  of  three  neurons. 

The  spontaneous  firings  of  neurons  A  and  B  form  independent  processes 
N^(*)  and  A  stimulus  neuron  (neuron  S)  has  spikes  which  form  the 

process  N^O).  Let  X^  equal  the  rate  of  N^(*),  for  i=l,2,3.  Suppose 
that  every  spike  from  S  deterministically  gives  rise  to  a  spike  from  A 
and  a  spike  from  B  t  and  £g  units  of  time  later,  where  and  Jig 

are  fixed  constants.  Figure  1  gives  a  diagram  depicting  this  situation. 

The  overall  spike  trains  of  neurons  A  is  therefore  the  superposition  of 
N^(*)  and  delayed  by  and  similarly  for  neuron  B.  We  there¬ 

fore  have 
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If  in  addition  N^(*)  is  a  Poisson  process,  a  simple  calculation  gives 


(7) 


XA  if  U  *  V£B 

\ 

if  U  =  £a"*b 


Figure  1:  Neuron  S  excites  A  and  B  with  latencies 
and  £g,  respectively. 

When  XA|g(*)  is  not  constant,  the  processes  NA(*)  and  Ng(*)  are 
dependent.  It  should  be  stressed  that  if  this  is  the  case,  we  cannot  infer 
causality:  neither  A  nor  B  need  be  affecting  the  other  directly.  Instead 
they  may  both  be  affected  by  a  third  source,  as  Example  2  illustrates. 


Since 


',!„(•)  can  indicate  deviations  from  independence,  it  is 

A  |  D 

important  to  be  able  to  estimate  it  from  the  data.  Cox  (1965),  Cox  and 

Lewis  (1972)  and  Brillinger  (1976)  have  proposed  estimates  of  the  function 

a a | B ( * ) »  or  what  is  essentially  equivalent,  Xg  •  Their  estimate 

is  formed  as  follows:  From  the  processes  <  •••  <  A  , 

A 

B.  <  B„  <  • • •  <  B„  ,  compute  all  the  differences  A.-B.,  for  1  <  i  <  N,, 
12  Ng  l  j  —  —  A 

i  <_  j  Then  form  a  histogram  (or  more  generally  any  density  estimate) 

of  these  N,  •  N’  points.  This  histogram  is  called  a  cross-correlation 
histogram  (CCH) . 

Brillinger  (1976)  showed  that  if  the  processes  (A^ }  and  {B^}  are 

observed  over  a  period  of  length  T,  then  as  T  -+■  00  very  roughly,  the 

suitably  normalized  CCH  resembles  the  function  '■!„(•).  In  a  little 

A  I  D 

more  detail,  suppose  that  the  bin  width  used  to  form  the  CCH  is  b.  His 

result  is  that  under  some  regularity  conditions,  if  T  ■*  «  and  b  ■+  0 

in  such  a  way  that  bT  remains  constant,  then  the  height  of  the  normalized 

CCH  at  a  fixed  point  u  has  mean  and  variance  that  is  of  the 

order  of  .  His  result  applies  to  individual  points  u,  and  does  not 
b  1 

imply  that  the  normalized  CCH  as  a  whole  resembles  ^Ib^*^’  furthermore, 
the  variance  of  the  height  of  the  normalized  CCH  is  of  a  larger  order  of 
magnitude  than  Nevertheless,  his  result  aids  greatly  in  understanding 

the  way  in  which  the  normalized  CCH  estimates 

We  now  illustrate  the  use  of  the  CCH's  in  the  situations  given  by  the 
two  examples  above. 

Example  A:  Two  independent  neurons. 

This  is  the  situation  described  in  Example  1.  The  parameters  are 
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and 


X  =  x2  =  100  , 

T  =  1 

(recall  that  T  is  the  length  of  the  period  of  observation).  Two  features 
of  the  CCH  are  apparent.  First,  as  expected,  the  CCH  is  generally  flat,  or 
at  least  there  are  no  apparent  peaks  or  valleys.  Second,  although  the 
heights  of  the  different  bins  in  the  CCH  all  have  the  same  distribution, 
they  clearly  display  a  lot  of  fluctuation,  in  accordance  with  the  remark 
made  concerning  the  variance. 

Example  B:  Two  neurons  simultaneously  excited  by  a  third  source. 

This  is  the  situation  described  in  Example  2.  The  parameters  are 

X1  =  X2  =  X3  =  100 

T  =  1 

eA  -  .04 

S  '  -03  • 

In  this  situation  we  have  (see  equation  (7)) 

fx.  if  u  +  .01 

A 

(S)  xA|B(u)  -  , 

»  if  u  =  .01  . 
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ntly . 


Figure  3  gives  the  CCH  for  this  setup.  The  location  of  the  spike  is  at 
u  =  .01,  in  accordance  with  (8). 

Example  C:  Two  neurons  simultaneously  excited  by  a  third  source. 

This  is  the  same  as  Example  B,  except  that  T  =  10.  The  function 

X.  !_.(•)  is  still  given  by  (8).  Figure  4  gives  the  CCH  for  this  situation. 
A  |  D 

The  spike  is  much  sharper  than  in  Figure  3,  due  to  the  fact  that  T  is 
much  larger. 
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150- 
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0001 
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neurons  simultaneously  excited  by  a  third  source 
f  observation  is  longer  than  in  Figure  3. 


3.  ESTIMATES  AND  CONFIDENCE  INTERVALS  FOR  THE  NORMALIZED  RENEWAL  FUNCTION. 

In  Example  2,  the  lags  £A  and  £_  were  taken  to  be  fixed  constants, 
and  this  led  to  the  highly  discontinuous  form  of  \^ | ^ ( • )  given  by  (8). 

We  saw  that  the  CCH  was  particularly  suited  for  revealing  the  dependence 
structure  in  this  situation.  In  practice  however,  the  lags  will  never  be 
fixed  constants;  this  will  result  in  being  a  smooth  function. 

The  normalized  CCH  will  then  not  be  well  suited  for  estimating  it:  it  is 
more  suitable  as  a  pointwise  rather  than  as  a  global  estimate  of 

Example  3:  A  two-neuron  network. 

Consider  the  network  described  in  Example  2,  except  that  =  0 
and  £g  is  random.  More  specifically,  if  denotes  the  lag  between 

impulse  i  of  the  stimulus  neuron  and  the  induced  impulse  from  neuron  B, 
assume  that  the  are  iid  with  a  density  f.  Then  it  is  simple  to 

D 

see  that 

<9>  SA|B(U)  ’  ■ 

In  this  example,  suppose  for  simplicity  that  f  was  the  standard 
uniform  distribution  on  [0,1].  Then,  if  a  CCH  was  formed  it  would  be 
likely  to  have  a  peak  between  0  and  1.  If  we  knew  the  distribution 
of  the  area  under  the  histogram,  we  would  be  able  to  determine  whether  the 
peak  was  significant  or  rather  was  due  to  random  fluctuation. 

Let  AR(tpt2)  denote  the  area  under  the  histogram  between  the 
points  t^  and  ^2 •  Here  t^  and  t2  are  arbitrary  points  satisfying 
t-^  <  t^ •  We  may  then  write 
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(10) 


AR(t1,t2)  =  b  l 
i=l 


I  I{Ai_Bi e  (t,,t2)}  , 

j=l  J  ' 


where  b  is  the  bin  width. 

Let  us  now  define  for  t^  <  t2 


(11) 


UA|B(tl’t2)  E{NA(t1»t2) 


There  is  a  B  point-, 
at  t  =  0  i 


UA | B ( ' » * )  is  called  the  conditional  renewal  function  or  renewal  function, 
for  short.  Note  that  if  N^(-)  and  Ng(*)  are  independent  processes, 
then  U^|  b^I1  =  ^A^t2-tl^  by  Korolyuk's  theorem.  Defining  the 
normalized  renewal  function  W(*,«)  by 


(12) 


W(tl’t2)  =  X  UA|B(tl’t2) 
A  ' 


* 


we  then  have 


(13) 


W(trt2)  = 


t2“tl 
>  t2-t]L 


<  t2“tl 


(independence) 

(excitation) 

(inhibition) 


with  the  words  "excitation"  and  "inhibition"  suitably  referring  to  the 
interval  (t^.t^).  Thus,  the  family  of  parameters  (W(t^,t2);  t^  <  t2> 
is  useful  in  describing  the  dependence  structure  between  N’A(»)  and 
Ng  ( * )  . 
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An  appropriate  way  to  search  for  an  excitatory  or  an  inhibitory 
effect  is  to  proceed  as  follows.  Fix  some  constant  A  >  0,  and  consider 
W(t  -  t  +  ^-)  as  a  function  of  t.  Here,  A  is  determined  by  the 
experimenter  as  the  likely  duration  of  an  interaction,  and  is  determined 
by  physiological  considerations.  Under  independence  this  function  is 
constantly  equal  to  A,  so  that  deviations  from  A  indicate  a  dependence 
structure.  It  is  therefore  necessary  to  estimate  W(t  -  -j,  t  +  -j) . 

In  Doss  (1983)  it  is  shown  that  under  some  regularity  conditions, 

AR(t  -  y,  t  + -|),  suitably  normalized,  is  asymptotically  normal,  with  mean 

A  A  1 

W(t  -  y,  t  +  -j)  ,  and  variance  of  the  order  of  More  specifically,  the 

result  is  that  for  large  Nfi, 

2 

(i)  AR(t  -  -|,  t  +  -|)  (b~NT~)  is  approximately  N(W(T  -  -|,  t  + -|) ,  °  -^) . 

BA  B 

Furthermore, 

2 

(ii)  o  (t)  can  be  estimated  consistently  from  the  data. 

The  main  feature  of  the  result  is  (ii),  which  allows  confidence  intervals  to 
be  formed.  This  makes  possible  a  formal  analysis.  The  entire  function 

AR(t  -  t  +  y)  -  can  be  plotted,  and  a  confidence  band  can  be  put 

A  B 

around  it.  Under  the  hypothesis  of  independence,  the  function  is  essentially 
flat,  at  height  A.  Upwards  or  downwards  deviations  from  A  (dependence) 
can  therefore  be  discerned  at  a  glance  from  the  plot. 
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Example  D:  Two  neurons  simultaneously  excited  by  a  third  source. 

This  is  the  situation  described  in  Example  2.  The  parameters  are: 

A1  =  X2  =  X3  =  100»  C  =  3  » 

=  .1,  =  .05,  (thus  £  is  negative) 

and  A  =  .1  . 

Figure  5  shows  the  estimate  of  the  normalized  renewal  function 

W(t  -  y,  t  +  -|)  for  -.5  <_  t  £  .5  . 

The  peak  in  the  diagram  clearly  shows  the  dependence  structure. 

In  Figure  5,  the  estimate  of  the  normalized  renewal  function  appears 
without  the  confidence  band.  Also,  the  diagram  was  constructed  from  data 
where  the  lags  were  fixed  and  not  random.  In  a  later  paper  we  hope  to  carry 
out  the  analysis  further  by 

(i)  constructing  the  confidence  bands 

(ii)  using  data  generated  as  in  Example  3,  i.e.,  with  random  lags,  and 

(iii)  illustrate  all  the  procedures  discussed  here  on  real  data  as  well. 
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for  the  second.  Such  data  arises  in  many  contexts,  and  it  is  often  of 


interest  to  discover  and  quantify  the  association  between  the  two  processes. 

Two  fields  in  which  this  situation  occurs  are  neurophysiology  and  reliability 
theory.  In  this  article,  we  describe  and  discuss  certain  graphs,  plots,  as 
well  as  more  formal  methods  that  can  assess  the  dependence  between  Doint 
processes.  Specifically,  these  methods  indicate  whether  or  not  the  likeli¬ 
hood  of  an  A-point  is  increased  (decreased)  after  the  occurrence  of  a 
B-point.  The  techniques  are  illustrated  on  simulated  data.  Although  bivariate 
point  processes  arise  in  many  fields,  we  emphasize  the  applications  in 
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